Crank-Nicolson demo

This implementation is again in Julia. Since symmetric tridiagonal matrices are so prevalent and important, Julia has a special data type for them. This makes implementing Crank--Nicolson simple. If you are using python or Matlab you want to use the spdiags command. For example, in python:

import numpy as np
from scipy.sparse import spdiags

m = 4
data = np.array([np.ones(m), -2.0*np.ones(m) , np.ones(m)]);
diags = np.array([-1, 0, 1])
A = spdiags(data, diags, m, m)

A.toarray()  # just to see what it looks like

Initial and boundary data

$$ \begin{cases} u_t = u_{xx},\\ u(x,0) = \eta(x),\\ u(0,t) = g_0(t),\\ u(1,t) = g_1(t). \end{cases} $$

Create a gif

Real-time plotting

Timing without plotting

Forward Euler-like discretization (what not to do...)